Dynamics of two-dimensional coherent structures in nonlocal nonlinear media 
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We study stability and dynamics of the single cylindrically symmetric solitary structures and 
dipolar solitonic molecules in spatially nonlocal media. The main properties of the solitons, vortex 
solitons, and dipolar solitons are investigated analytically and numerically. The vortices and higher- 
order solitons show the transverse symmetry-breaking azimuthal instability below some critical 
power. We find the threshold of the vortex soliton stabilization using the linear stability analysis 
and direct numerical simulations. The higher-order solitons, which have a central peak and one or 
more surrounding rings, are also demonstrated to be stabilized in nonlocal nonlinear media. Using 
direct numerical simulations, we find a class of radially asymmetric, dipole-like solitons and show 
that, at sufficiently high power, these structures are stable. 

PACS numbers: 42.65.Sf, 42.65.Tg, 42.70.Df, 52.38.Hb 



I. INTRODUCTION 

The recent experimental observations of spatial soli- 
tons in nonlocal media such as nematic liquid crystals 
PJ , lead glasses @ , renewed an interest to coherent struc- 
tures in spatially nonlocal nonlinear media. In the spa- 
tially nonlocal media the nonlinear response depends on 
the wave packet intensity at some extensive spatial do- 
main. Nonlocality is a key feature of many nonlinear 
media. It naturally appears in different physical systems 
such as plasmas |1J5[ , Bose-Einstein condensates (BEC) 
, optical media Q , atomic nuclei @ , liquid crystals Q . 

An important property of spatially nonlocal nonlinear 
response is that it prevents a catastrophic collapse which 
usually occurs in local sclf-focusing media when the 
power of the two- or three-dimensional wave packet ex- 
ceeds some critical value. In particular, a rigorous proof 
of absence of collapse during the wave packet propagation 
described by the nonlocal nonlinear Schrodinger equation 
(NLSE) with sufficiently general symmetric real-valued 
response kernel has been presented in 0, EJ ■ In the ab- 
sence of collapse, the competition between diffraction 
spreading and nonlinear self-action leads to formation 
of the stationary solitary wave structures - solitons and 
vortex solitons. Different types of two-dimensional self- 
trapped localized wave beams have been predicted and 
experimentally observed in various nonlinear media |lfjj |. 

Fundamental solitons are the lowest-order localized 
structures with a single peak in the intensity distribu- 
tion and transversely constant phase. A vortex is the 
structure with ring-like intensity distribution, with the 
dark hole at the center where the phase dislocation takes 
place: a phase circulation around the axis of propagation 
is equal to 2irm. An integer m is refered to as topologi- 
cal charge. The important integral of motion associated 
with this type of solitary wave is the angular momen- 
tum, which can be expressed through the vortex power 
N and topological charge m as follows: |M| = mN. 
Thus, vortices (or spinning solitons) are the localized 
structures with nonzero angular momentum. While the 



fundamental solitons are robust in a collapse-free me- 
dia, the spinning solitons may possess a strong azimuthal 
modulational instability As a result, the vortex de- 
cays into several ordinary solitons which fly off carrying 
out its energy and angular momentum. Though vortices 
can be stabilized in media with competing local nonlin- 
caritics 

[3 111 111 113 , it occurs only in the extreme 
regime when the higher-order contribution to the nonlin- 
earity dominates. The recent investigations [l6l[T7| have 
shown that the different kinds of nonlocality of the non- 
linear response can also suppress or completely eliminate 
the symmetry-breaking azimuthal instability for one-ring 
and double-ring single-charge (m = 1) vortex solitons. 
The recent experiments Q have confirmed an existence 
of stable single-charge vortex solitons in the the nonlo- 
cal media with thermal optical nonlinearity. The multi- 
charge vortices (m > 1) are unstable in the media with 
thermal nonlocal nonlinearity, as was shown in Ref. 0] 
by linear stability analysis of small azimuthal perturba- 
tions and direct numerical modeling. On the other hand, 
the robust propagation of two-charge vortex (m = 2) 
solitons has been observed in the numerical simulations 
[l7j for the model based on Gaussian-type kernel of the 
nonlocal media response function. We perform here the 
linear stability analysis for the model used in Ref. 0] 
and prove the possibility of stabilization of multicharge 
vortex solitons in the nonlocal nonlinear media. 

The higher-bound solitons with field nodes (zero- 
crossing) have been first discovered in Ref. [2(| for the 
local Kerr-type nonlinear media. The nth bound state 
has a central bright spot surrounded by n rings of varying 
size. In the local nonlinear media the higher-order soli- 
tons with zero angular momentum show the azimuthal 
instability [2?l l2Sj similar to the instability of the vor- 
tex solitons. The rings which surround the central peak 
possess an azimuthal instability. As a result, the higher- 
bound structures decay into several fundamental solitons. 
We reveal here that nonspinning higher-order solitons can 
be stabilized in the nonlocal medium. 

Another important feature of nonlocal nonlinear me- 
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dia is the possibility of existence of composite soliton 
structures. A composite soliton structure, or a multi- 
soliton complex is a self-localized state which is a nonlin- 
ear superposition of several fundamental solitons 0, . 
Multi-soliton structures in nonlocal media were consid- 
ered first in Refs. 0, |2{1, and they have recently re- 
ceived renewed interest |2ll l22l l23j . In particular, one- 
dimensional (ID) nonlocal model suggested in Ref. Q| 
was studied in Ref. and it was shown that dipole-, 
triple-, and quadrupole-mode solitons can be made sta- 
ble. Very recently [24| two-dimensional (2D) rotating 
dipole structures were considered in the framework of an 
approximate variational approach. In this paper, using 
direct 2D simulations, we find numerically a class of ra- 
dially asymmetric two-dimensional dipole-mode soliton 
solutions and show that, at sufficiently high input power, 
these solutions are stable. 

The aim of this paper is to study general properties 
and to carry out the stability analysis of single solitons 
(both fundamental and higher-bound states) , vortex soli- 
tons, and composite dipole-like solitons in the strongly 
nonlocal media. 

The paper is organized as follows. In Sec. ^ we formu- 
late a model and present basic equations. The stability 
analysis based on the variational approach and numerical 
simulations for single solitary structures is performed in 
Sec. IIIII and then, in Sec. IIVI we consider the multisoli- 
ton dipole-like structures. The conclusions are given in 
Sec. E| 



The nonlocal nonlinear media response function can be 
taken as follows: 

0(f) = J Ri^-riD^irDfdW (3) 

The shape of the kernel R(r) is determined by the type 
of the nonlocal interaction in media and can be rather 
complicated Il6j. However, there are general proper- 
ties valid for all nonlocal media response functions. The 
nonlinear term tends to the local Kcrr-type nonlincarity: 
— ► | "J | 2 when the spatial scale of the wave packet in- 
tensity distribution \^>\ 2 is much wider than the effective 
width of the potential R(r). In the opposite case of the 
strongly nonlocal regime, the response function can be 
estimated as follows: 0(r) = N {R(0) + ±Ai?(0)r 2 }. In 
the latter case, the highly-nonlocal NLSE (1} with a suffi- 
ciently regular kernel R(r) is mathematically identical to 
the linear Schrodingcr equation with harmonic oscillator 
potential, as was pointed out in Ref. [29j . 

We consider in this paper the nonlocal response func- 
tion kernel modeled by the Gaussian shape kernel: 

2 

R Qf-f[\) = SL e-aV-rlf (4) 
7T 

where a is the nonlocality parameter. Keeping the main 
features of nonlocal media this model allows a very ac- 
curate and simple analytical treatment. 

III. SINGLE SOLITARY STRUCTURES 



II. BASIC EQUATIONS 

We consider propagation of the electric-field envelope 
^/(x, y, z) described by the paraxial wave equation: 



9* 

i— + + = 0, 

uz 



(1) 



where z is the direction of beam propagation, repre- 
sents the nonlinear response of the media. 

Equation £Q) conserves the following integrals of mo- 
tion: (i) number of quanta ("energy" or "beam power"): 



N = I m 2 d 2 r, 



(2) 



(ii) momentum: 

I_L = ~ J (**V_l* - *V_l**) d 2 r, 

(iii) angular momentum: 

M = ~ J [r x (**V ± * - tfVj.**)] d 2 r, 

(iv) Hamiltonian: 



|V^*| 2 --0|*| 2 J.d 2 r. 



In this section, we study single solitons (both funda- 
mental and higher-bound states) and vortex solitons. We 
look for the stationary solutions of the Eq. (JTJ in the 
form: 



$(x,y,z) = i(j(r)e 



imip-\-iAz 



(5) 



where (p and r = \J x 2 + y 2 are the azimuthal angle and 
the radial coordinate, respectively, and A is the beam 
propagation constant. Such solutions describe cither the 
soliton, when m = 0, or the vortex soliton with the topo- 
logical charge m, when m ^ 0. The function ip(r) obeys 
the integro-differential equation: 



(0) 



The boundary conditions for the localized solutions are: 
ip(oo) = and (dtp/dr) r=0 = for solitons, -0(0) = 
for vortices. For the stationary solution JSJ it is easy to 
rewrite the response function in the form: 

r+oo 

6{r) = 2a 2 / e- a2{r - ri)2 l (2a 2 rr 1 )\iP(r 1 )\ 2 r 1 dr 1 , 
Jo 

(7) 

where X v (x) = e x I v (x) is the exponentially scaled mod- 
ified Bessel function. 

In the next subsection we start our considerations on 
the single cylindrically symmetric solitary structures with 
the analytical analysis based on the variational approach. 
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A. Variational approach 

As known j2|j, the nonlocal NLSE turns to the linear 
Schrodinger equation with harmonic oscillator potential 
in the highly-nonlocal limit, when the spatial scale of the 
response function is much wider than the wave packet 
localization region. Since the Laguerre-Gauss modes are 
the exact eigenstates for the two-dimensional linear os- 
cillator, the variational method with the trial function of 
the form 

*(x,y,z) = h(z)C n L^ (£ 2 ) e -ie{i+Mz)}+*rn V+l $(z)^ 

. (8) 

is expected to give an accurate description of all eigen- 
modes, especially in the highly nonlocal regime. Here 
L^ n \x) is the generalized Laguerre polynomial, n is 
the number of nodes of the radial profile, and m is 
the topological charge, £ = r/a(z), where a(z) is the 
first variational parameter that characterizes a radius 
of solitary structure: a 2 = (r 2 )(2n + m + where 
(r 2 ) = A _1 (^ , |r 2 |\E') is the mean-square radius. The sec- 
ond variational parameter b(z) is the phase curvature. 
The amplitude h(z) can be readily found from the rela- 
tion h(z)a(z) = \J n(n+iitf i ' wrucn is obtained from nor- 
malization condition |2f . 

We start our considerations with the lower-order node- 
less states (n = 0). The nonlinear response function in 
the framework of variational approach is given by expres- 
sion: 

0(r) - h\e-^m\{l - q) m L m ( q 2 e/(q - 1)) , 

where q = a 2 a 2 /(l + a 2 a 2 ), L m (x) = L„ (i) is Laguerre 
polynomial of m-th order. Note, if a 2 a 2 ^S> 1, then q — > 1, 
and one obtains — > \^\ 2 , as it should be in the local 
limit. 

In accordance with the variational procedure, we con- 
struct the Lagrangian density 



2 1 dz 



and the Lagrangian 



9* 

dz 



|v^| 2 _I e |vi/| 2 



C = I l d 2 r = A$ + (m+1) AT 



where b = b/a, H is the Hamiltonian: 
H = N I (m + 1) (l/a 2 + b 2 ) - 



ba — ba 



Nhmja) 
2ira 2 



where h m (a) = q{l- q) m (1 + q) m 1 P m (jz^i 
P m {x) is the Legendre polynomial of the m-th order. 

The first two dynamical equations can be written in 
the canonical form: 



N(m + 1) da 
2 dz 



dH 
~db' 



N(m + 1) db 
2 dz 



dH 
da 



(9) 



and the third one is the following: dN/dz = which 
means that the number of quanta is the integral of mo- 
tion. The soliton or vortex soliton corresponds to the sta- 
tionary point of the Hamiltonian: dH/db = 0, dH/da = 
0. The first condition yields &o = 0. From the second 
equation one can find the width ao as the function of the 
number of quanta N . It is easy to verify that vortex or 
soliton can exist only above the threshold value for num- 
ber of quanta: N > N cr = 4 m+1 7r(m + l)(m!) 2 /(2m)!. 
Using the similar procedure we have considered the non- 
spining (m = 0) higher-bound states with one (n = 1) 
and two nodes {n = 2). The thresholds for an existence 
of the higher bound states are as follows: N cr = 24ir for 
n = 1 and N cr = 640tt/11 for n = 2. 

The described variational procedure provides the pos- 
sibility for analysis of the stationary radially-symmetric 
coherent structures. The results of the variational anal- 
ysis are given in Fig. ^ (a) for fundamental solitons and 
vortex solitons, and in Fig. [2 for higher-order solitons. 

Moreover, using the set © it is possible to study the 
radially-symmetric dynamics of the localized wavepack- 
ets propagating in z-dircction. Let us investigate, for 
example, the evolution of a slightly perturbed stationary 
soliton solution. It can easily be shown that the small 
deviation of the soliton width S = a — ao from the sta- 



tionary value ao obeys the equation 5 



0, where 



8N~ 1 (d 2 H/da 2 ) a0i b - Therefore, the soliton being 
radially perturbed exhibits the osci llations with the fre- 
quency u = 8a 2 7r^ 1 / 2 [l - ^Ncr/N} 3 / 2 , where N cr = 4tt 
is the threshold for soliton existence. 

However, one cannot study the stability of stationary 
solutions with respect to symmetry-breaking azimuthal 
perturbations in the framework of a variational approach 
with a radially symmetric trial function. Stability con- 
ditions of steady-state solutions, regarding small general 
2D perturbations, will be obtained by a linear stability 
analysis in the next subsection. 



B. Numerical modeling 

The boundary problem © is equivalent to the integral 
equation: 



«(#(f/)G ro (» l ,r;VA)r)(i ) ), (10) 



where 9 is given by Eq. J7J), 



Gm^i^a) 



4(a6)4«i), < & < 6, 
I m {a£2)K m (a£i), £ 2 < £l < +°o, 



(11) 

where I m and K m are the modified Bcsscl functions of the 
first and second kind, respectively. We have solved the 
nonlinear integral equation l|10|) using stabilized iterative 
method |25|. For numerical modeling it is useful to per- 
form the rescaling transformation of the form: r' = ar, 
ip' = tjj/a, 6' = 9 /a 2 , z' = za 2 . Such transformation 
reduces the number of parameters to one dimcnsionlcss 
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FIG. 1: (a) Number of quanta N vs parameter A = A/a 2 for solitons (n = 0,m = 0) and vortex solitons (n = 0, m = 1,2) 
(solid curves for variational and circles for numerical results). Numerically found profiles for (b) A = 0.05; (c) A = 0.5; (d) 
A = 5. Solid curves for ip'{r') = tp/a dashed curves for 0'(r') = 9/a 2 ; the scaled coordinates r' = ar are used. 



parameter A = A/a 2 . Figure ^ (b)-(d) shows several ex- 
amples of the numerical solution of the © at different 
values of the parameter A. Note, that the nonlocal limit 
a 2 <C A corresponds to the large values of the parameter 
A and, as seen from Fig. n( a )> to the large values of the 
beam power N. 

Let us investigate the stability of the steady-states with 
respect to small azimuthal perturbations. Expanding the 
nonstationary solution in vicinity of a steady-state: 



iAz-\-imtp 



#(r,z) = {ip(r) +Sip+ +<5?Ai} e 



where 5ip±(r,ip,z) = e±(r)e tu)z+tLv , and linearizing the 
dynamical Eq. one can obtain the eigenvalue problem 
for u> of the form: 



Qrn+L + 9L 

-9l 



9L 

Qrn-L — 9L 



e = use, 



(12) 



where e = (e+,e_), and |Imct>| determines the growth 
rate of an unstable mode, 

Qm± L e± = j— A + A( m±i > + 6(r)} e±, 



g L e± = 2a 2 i&(r) / i/>(t)x{r, 0e± 



x(r,£) = £,e~ a2< - r ~^ 2 lL(2a 2 r£ i ). The unperturbed radial 
profile ip(r) is taken to be real without loss of generality. 
We have used the Hankel spectral transformation and re- 
duced the integro-differential eigenvalue problem (|12ll to 
the linear algebraic one. Figure|21(a) shows the maximum 
growth rate for one-charge (m = 1) vortex solitons. It 
is seen that only modes having L = 1,2,3 can be unsta- 
ble. Modulational instability is strongly suppressed in a 
highly-nonlocal regime: the growth rates vanish at some 
finite values of parameter A. The mode with azimuthal 



number L = 2 corresponds to the largest growth rates 
with widest instability region: all growth rates are equal 
to zero at A > X cr sa 9.1. Similar analysis has been per- 
formed for multi-charge vortices (with m = 2, ...5). Fig- 
ure[s2(b) depicts the growth rates for two-charge (m = 2) 
vortex solitons. Note that the modulation instability is 
eliminated for A > X t h sa 23.8. Thus, in the media with 
the nonlocal response of the form the multi-charge 
vortices can be stabilized. Note that vortex solitons with 
m > 1 are found [l6l | to be unstable in the media with 
thermal nonlocal nonlinearity. Therefore, the dynamical 
properties of the vortex structures can be significantly 
affected by specific type of the nonlocal interaction. 

The results of linear analysis have been confirmed by 
extensive series of numerical simulations of dynamics of 
perturbed stationary solutions. We have performed the 
split-step Fourier transform method to solve the dynami- 
cal Eq. Q with the response function Eq. J3J). The non- 
local nonlinear term has been calculated in the spectral 
domain, since it has the form of the convolution of the in- 
tensity distribution ^(r)! 2 with the function i?(|r|). We 
have used the numerically found stationary vortex soli- 
tons and variational profiles for the higher-order solitons 
as the initial conditions. Different kinds of the perturba- 
tions such as random noise, radially-symmetric and az- 
imuthaly periodical perturbations have been applied in 
the numerical experiments. The conclusions of the linear 
stability analysis are found to be in a good agreement 
with our simulations: the vortex solitons become stable 
above some critical power which is close to the one pre- 
dicted by linear stability analysis. 

Figure 0] (a) illustrates the unstable evolution of the 
single-charge vortex soliton. In the stable region, when 
A > A cr , vortices survive even being strongly perturbed. 
The mean-square radius and intensity of the vortex os- 
cillate, but the vortex ring shows the robust propagation 
over vast distances (thousands of diffraction lengths) for 
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FIG. 2: (a) Number of quanta iV vs parameter A = A/a 2 for higher-order solitons (m = 0, n = 1,2). Radial profiles for (b) 
A = 0.05; (c) A = 0.5; (d) A = 5. Solid curves for ip' dashed curves for 8' . 




FIG. 3: (Color online) The scaled growth rate 7 = Imui/a 2 of 
linear perturbation modes vs the parameter A = A/ a 2 for the 
vortices with (a) m = 1 and (b) m = 2, the inset depicts the 
growth rates of the high-L modes in more details. The num- 
bers near the curves stand for the azimuthal mode numbers 
L. 



the hundreds of the effective periods of the oscillation. 
The higher-order solitons exhibit similar decay in weakly 
nonlocal regime. However, when beam power increases, 
an interesting dynamics with revivals has been observed 
[see Fig. H|(b)]. Initially field envelope decays into sev- 
eral filaments, but then it recurs at larger propagation 
distances. With further increasing of the power, a higher- 
order soliton occurs to be stabilized: it shows the robust 
propagation without breakup during the hundreds of ef- 
fective oscillation periods. 



IV. STABLE BISOLITON MOLECULES 

In this section, we present and study localized asym- 
metric dipole-like solutions. We look for stationary solu- 
tions of Eq. Q in the form \I/(cc, y, z) = ip(x, y) exp(iAz), 



z' =0 


z' =2 


z' 







• 




• 











(a) 





-1 1 


-1 1 


-1 1 




X' 


X' 


X' 




z' =0 


z' =1.8 


z' =3 


0.2 




• 




y 

-0.2 


© 


Y 






(b) 



-0.2 0.2 
x' 



-0.2 0.2 

X' 



-0.2 

X' 



FIG. 4: Snapshots of \V(x',y')\ for different z' . (a) Example 
of unstable evolution of single-charge (m = 1, n = 0) vortex 
soliton with A = 4. (b) Dynamics with revival of the non- 
spinning (m = 0) one-node (n = 1) soliton at A = 35. 



so that i/j(x, y) obeys the equation 

-Aip + A ± <ip + 9Tp = 0, 

where 



(13) 



£L / e -^ [{x - x i) a +(v-»0V ( )d a ri> (14) 

7T ' 



and we do not assume the radial symmetry of ip(x,y). 
Unlike the NLSE with local cubic nonlincarity, equation 
(|13f) with the nonlocal nonlinear media response 8 given 
by Eq. I|14|) has two characteristic transverse scales: the 
internal scale A 1 / 2 (in the NLS this scale determines the 
characteristic size of the soliton) and the " external" scale 
a which is the measure of nonlocality. Under this, the 
characteristic size of the self-consistent potential well in 
Eq. 0} can significantly differ from A -1 / 2 and, thus, the 
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FIG. 5: The dipole solution ip' = ip/a for A = A/a 2 = 40; 
the scaled coordinates x' = ax and y' = ay are used. 



existence of composite soliton structures becomes pos- 
sible. In this paper we restrict ourselves to dipole-likc 
localized solutions. Results concerning the tripolar and 
higher radially asymmetric soliton modes will be pub- 
lished elsewhere. 

As above, we use the scaled variables tp' , z', x' = ax, 
and y' = ay. Imposing periodic boundary conditions on 
Cartesian grid and choosing an appropriate initial guess, 
one can find numerically dipole-type localized solutions 
by using the relaxation technique similar to one described 
in Ref. j^. An example of such dipole solution is pre- 
sented in Fig. The dipole consists of two out-of-phasc 
monopoles. The characteristic width of the monopoles 
in the dipole and the "distance" between them decrease 
with increasing the parameter A = A/a 2 . 

We next addressed the stability of these dipole solu- 
tions and study the evolution (propagation) of the dipoles 
in the presence of small initial perturbations. We have 
undertaken extensive numerical modeling of Eqs. lfT|). Q 
and initialized with our computed dipole-type solu- 
tions with added gaussian noise. Spatial discretization 
was based on the pseudospectral method and "tempo- 
ral" z-discretization included the split-step scheme. The 
numerical simulations clearly show that the dipoles with 
A > Xth, where X t h ~ 21 is the threshold value, are sta- 
ble with respect to small initial noisy perturbations up to 
the maximum propagation distances used (of the order 
of z' = 3000). The stable propagation of the dipole is 
illustrated in Figs.[5{a), (b). Additionally, the stable dy- 
namics was monitored by plotting the z' dependencies of 
the averaged intensity J |^| 4 (ir/7V and mean-square ra- 
dius / r 2 \ip\ 2 dr/N. For stable propagation, these quan- 
tities undergo small oscillations near the equilibrium val- 
ues. Note, that dipoles with sufficiently large (com- 
pared to X t h) values of A survive over huge distances 
(many thousands of diffraction lengths) in the presence 
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X=22 » 1 z'=2500 








FIG. 6: Snapshots of \ip'(x' , y') \ for dipole propagation in the 
presence of small random perturbation for different A = A/a 2 , 
shown in the scaled (x' = ax, y' — ay) plane, at different 
moments z: (a), (b) - stable propagation; (c), (d) - destruction 
of the original dipole; (e), (f) - decay into two solitons. 



of quite significant perturbations. We performed a se- 
ries of runs for A > 200 in the presence of strong ini- 
tial noise. The initial condition was taken in the form 
ip'(x', y')[l+ef(x', y')], where ijj'(x',y') is the numerically 
calculated exact dipole solution, f(x',y') is the white 
gaussian noise with variance a 2 — 1 and the parameter 
of perturbation e = 0.1 -j- 0.3. Snapshots of (x' ,y')\ 
at different z' for the case A = 400 and e = 0.12 are 
presented in Fig. [7| One can see that the dipole turns 
out to be extremely robust - even at z' = 2000 one can 
not detect any substantial distortion of the dipole shape. 
The dipoles, however, become unstable (even if the ini- 
tial noise is very small) if A < X t h- The typical decay 
of the unstable dipole near the threshold value of the 
rescaled propagation constant is shown in Figs.|ljc), (d). 

The situation, however, changes below X cr f=a 7.6. Un- 
der this, the dipole splits in two monopoles which move 
in the opposite directions without changing their shape, 
i. e. the monopoles just go away at infinity. In Fig. [S] 
we plot the dipole energy N<n p and the doubled energy 
2-ZV mon , where N mon is the energy of the monopole soli- 
ton solution of Eq. I|13f) . calculated numerically, versus 
the propagation constant A. One can see that the bound 
energy SN = Nd ip — 2N mon in the dipole tends to almost 
zero as A approaches X cr = 7.6. This explains why the 
dipole with A < A cr can be easily (i. e. under the action 
of extremely small initial perturbations) split into two 
monopole- type solitons. 

The results of the numerical simulation can be illus- 



FIG. 7: (Color online) Evolution of the dipole solution with 
A = 400 in the presence of strong initial random perturbation; 
the scaled variables x' = ax,y' = ay and z = o?z are used. 



FIG. 9: Contour maps of the function £(/3, fx) : (a) A = 40; 
(b) A = 7.6 
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FIG. 8: The doubled energy of the soliton (monopole) 2N mon 
(solid line) and the energy of the dipole Ndi P (dashed line) v. s. 
the parameter A = A/a 2 . 



trated through the variational analysis. Equation 1131) 
(in the scaled variables) is the Euler-Lagrangc equation 
for the Lagrangian 



£ = 



|W'| 2 - ^ ,2 + A^' 2 



d 2 r. (15) 



Taking a trial function in the form 



y/ = Ae -?[{*'-d/v 2 +v'*\ _ Ae -f-[( X '+d/2fw 2 ]^ (16) 



where A, f3 and d are unknown parameters to be deter- 
mined by the variational procedure, and substituting it 
into Eq. lfTo|l . we get 



C = nA z [ 2ci - ^-c 2 + Ac 3 I , 



(17) 



where 



ci = l-(l-M 2 /2)e^ / 2 , 
1 



C2 



■[l + 2e 



-M 



/? 2 (1 + /? 2 ) 

_ 4e -(3+2/3 2 ) M 2 /4(l+/3 2 ) 



(18) 

(19) 
(20) 



and, instead of d, we have introduced the variational pa- 
rameter [i = (3d. The optimum A satisfies the equation 
dC/dA = which yields 



A v 



4c x + 2Ac 3 

C2 



(21) 



The Lagrangian Eq. (|17fl . where A 2 is defined by Eq. 
I|2ip. depends only on two unknown variational parame- 
ters (3 and [i (or, cquivalently, f3 and d) and can be easily 
analyzed numerically. The topography of the function 
C(P, /x) depends on the rescaled propagation constant 
A = A/a 2 . There is the only minimum if A > A cr , where 
\ cr is some critical value. The contours (level lines) of 
the function £(/3, /x) for A = 40 arc shown (in the vicinity 
of the minimum) in Fig. [5] (a). In this case the mini- 
mum takes place at f3 — 1.74 and /x = 1.1. and corre- 
sponds to the dipole solution presented in Fig. |S| The 
amplitude of the approximate analytical solution A cal- 
culated from Eq. (|21|l and the parameters (3 and /x are in 
agreement with corresponding values estimated from the 
exact numerical solution - the comparison of the vari- 
ational analysis and the direct numerical simulation is 
presented in Fig. ^| (a). The dependence of the width 
of the monopoles in the dipole and the "distance" d 
between them on A is shown in Fig.^H(b). The topog- 
raphy of the function C{(3, fi) in the vicinity of the min- 
imum represents a long narrow valley oriented at some 
angle to the /x-axis. Under this, the depth of the valley 
and the angle to the /x-axis decreases with decreasing A. 
A similar situation holds for all A > A cr and we found 
A cr w 7.6. The picture changes sharply at A ~ X cr . The 
local minimum disappears, the saddle point arises, and 
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20 40 60 80 



FIG. 10: (a) The dependence of the dipole amplitude A on 
A): solid curve - variational analysis; dashed curve - direct 
numerical simulation, (b) Dipole distance d and the width 



s. A (variational analysis) 



this corresponds to the unstable dipole decaying into two 
monopolcs. The contours of the function C((3,fi) in this 
case (for A = 7.6) are shown in Fig. El (b). The found 
critical value X cr = 7.6 is in perfect agreement with the 
results of direct numerical simulation (see Fig.|Sfb)). 



V. SUMMARY AND CONCLUSIONS 

We have investigated the main properties and stability 
of the stationary two-dimensional localized solitary struc- 
tures in the nonlocal nonlinear media. We have studied 
both fundamental and higher-order solitons; one-charge 
and multi-charge vortex solitons with nonzero angular 
momentum; dipolar multisolitons. While the fundamen- 



tal soliton is always stable, the vortex solitons possess 
a strong azimuthal instability which is eliminated only 
in the strongly nonlocal regime. We have performed 
the linear stability analysis and direct numerical simu- 
lations to investigate the stability of vortices with arbi- 
trary topological charge. We have found the edge of the 
modulational instability and predicted the threshold for 
the beam power of the robust vortex soliton. We prove 
that in contrast to the nonlocal media with thermal non- 
linearity, the nonlinear response with the Gaussian-type 
kernel can sustain not only single-charge but also multi- 
charge vortices. We have investigated nonlocal higher- 
order nonspinning solitons which are the structures with 
the intensity distribution in the form of a bright spot sur- 
rounded by the bright rings. We theoretically predict an 
existence of stabilized higher-order nonspining solitons 
in the nonlocal media. Finally, we have found station- 
ary dipole-likc multisolitons which are the bound states 
of the out-of-phase fundamental solitons. We have simu- 
lated numerically the dynamics of the multisolitons in the 
presence of initial noise and performed simple variational 
analysis. It turns out that multisolitons are extremely ro- 
bust at sufficiently high input power in a highly nonlocal 
media. Therefore, these predictions open the prospects 
for the experimental observations of a wide class of stable 
coherent structures in various nonlocal nonlinear media. 
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